# Clearing environtment
rm(list = ls())

# Loading Packages
library(tidyverse)
library(magrittr)

# Loading Data
dat <- read.csv(file = "./RILA_Data_final.csv")

# Save options
sa <- 1
width <- 32
height <- 16

width2 <- 36
t1 <- 45
t2 <- 45

# Overleaf Path
opath <- "../../../Apps/Overleaf/Target Date RILAs/"

# Cleaning Data -----------------------------------------------------------
# Removing too short
dat %<>%
    filter((finish - start) / 60 > 30) %>%
    filter(grepl("385", la10) &
        grepl("200", la11))

# Retirement account only
dat <- dat %>% filter(aq6 == "Yes")

# l <- which(table(dat$Worker_ID)>1)
# k <- which(dat$Worker_ID %in% names(table(dat$Worker_ID))[l])
#
# dat <- dat[-k,]


# Updated function to calculate confidence interval
calc_ci <- function(p, n) {
    se <- sqrt(p * (1 - p) / n)
    ci <- 1.645 * se
    return(ci)
}

  # Financial Literacy ------------------------------------------------------
  dat$fin_lit <- if_else(dat$aq8 %in% c("Somewhat sure â€“ I understand most of what I will need to know", "Very sure â€“ I understand money management very well"),
                         1, 0)
  
  dat$fin_lit_savings <- if_else(dat$aq15 == "More than $102", 1, 0)

  dat$fin_lit_inflation <- if_else(dat$aq16 == "Less than today", 1, 0)

  dat$fin_lit <- dat$fin_lit_savings * dat$fin_lit_inflation

  summary(dat$fin_lit)

cor(dat$group_assigned, dat$gross_expense_ratio_random)
dat$change1 <- ifelse(dat$q1_a == dat$q1_b, 0, 1)
dat$change2 <- ifelse(dat$q2_a == dat$q2_b, 0, 1)
dat$change3 <- ifelse(dat$q3_a == dat$q3_b, 0, 1)
dat$change4 <- ifelse(dat$q4_a == dat$q4_b, 0, 1)

# Labeling Answers
dat$q1_a <- case_when(
    dat$q1_a == "Fund A" ~ "Mid Risk",
    dat$q1_a == "Fund B" ~ "More Risk",
    dat$q1_a == "Fund C" ~ "Less Risk"
)

dat$q2_a <- case_when(
    dat$q2_a == "Fund A" ~ "Mid Risk",
    dat$q2_a == "Fund B" ~ "More Risk",
    dat$q2_a == "Fund C" ~ "Less Risk"
)

dat$q3_a <- case_when(
    dat$q3_a == "Fund A" ~ "Mid Risk",
    dat$q3_a == "Fund B" ~ "More Risk",
    dat$q3_a == "Fund C" ~ "Less Risk"
)

dat$q4_a <- case_when(
    dat$q4_a == "Fund A" ~ "TDF",
    dat$q4_a == "Fund B" ~ "Floor",
    dat$q4_a == "Fund C" ~ "Buffer"
)


dat$group_assigned <- case_when(
    dat$group_assigned == 1 ~ "1. Base",
    dat$group_assigned == 2 ~ "2. Default Option",
    dat$group_assigned == 3 ~ "3. Distributional Info",
    dat$group_assigned == 4 ~ "4. Dist. and Default"
)

(log(.8) - log(.4)) * .12



btab <- dat %>%
    mutate(
        Male = as.numeric(aq1 == "Male"),
        White = as.numeric(aq2 == "Caucasian – not Hispanic" |
            aq2 == "Caucasian â€“ not Hispanic"),
        Married = as.numeric(aq3 == "Married"),
        HS = as.numeric(aq4 == "High school graduate or GED equivalent")
    ) %>%
    group_by(group_assigned) %>%
    summarize(
        n = n(),
        age_current_mean = mean(age_current),
        age_current_sd = sd(age_current),
        salary_mean = mean(salary_current / 1000),
        salary_sd = sd(salary_current / 1000),
        Male_mean = mean(Male) * 100,
        Male_sd = sd(Male) * 100,
        White_mean = mean(White) * 100,
        White_sd = sd(White) * 100,
        Married_mean = mean(Married) * 100,
        Married_sd = sd(Married) * 100,
        HS_mean = mean(HS) * 100,
        HS_sd = sd(HS) * 100
    )



# Balance Table
mdif <- function(m1, m2, s1, s2, n1, n2) {
    (m1 - m2) / (sqrt(s1 / n1 + s2 / n2))
}

# Formatting
f <- function(x, digits = 2) {
    if (x > 100) {
        format(x, digits = digits, nsmall = 0, big.mark = ",")
    } else {
        format(x, digits = digits, nsmall = 2, big.mark = ",")
    }
}

cat("Current Age & ",
    f(btab$age_current_mean[1]), "&",
    f(btab$age_current_mean[2]), "&",
    f(btab$age_current_mean[3]), "&",
    f(btab$age_current_mean[4]), "\\\\ & (",
    f(btab$age_current_sd[1]), ") & (",
    f(btab$age_current_sd[2]), ") & (",
    f(btab$age_current_sd[3]), ") & (",
    f(btab$age_current_sd[4]), ") \\\\ \\\\
      Average Salary (\\$000s) & ",
    f(btab$salary_mean[1]), "&",
    f(btab$salary_mean[2]), "&",
    f(btab$salary_mean[3]), "&",
    f(btab$salary_mean[4]), "\\\\ & (",
    f(btab$salary_sd[1]), ") & (",
    f(btab$salary_sd[2]), ") & (",
    f(btab$salary_sd[3]), ") & (",
    f(btab$salary_sd[4]), ") \\\\ \\\\
      Percent Male & ",
    f(btab$Male_mean[1]), "&",
    f(btab$Male_mean[2]), "&",
    f(btab$Male_mean[3]), "&",
    f(btab$Male_mean[4]), "\\\\ & (",
    f(btab$Male_sd[1]), ") & (",
    f(btab$Male_sd[2]), ") & (",
    f(btab$Male_sd[3]), ") & (",
    f(btab$Male_sd[4]), ") \\\\ \\\\
      Percent White & ",
    f(btab$White_mean[1]), "&",
    f(btab$White_mean[2]), "&",
    f(btab$White_mean[3]), "&",
    f(btab$White_mean[4]), "\\\\ & (",
    f(btab$White_sd[1]), ") & (",
    f(btab$White_sd[2]), ") & (",
    f(btab$White_sd[3]), ") & (",
    f(btab$White_sd[4]), ") \\\\ \\\\
      Percent Married & ",
    f(btab$Married_mean[1]), "&",
    f(btab$Married_mean[2]), "&",
    f(btab$Married_mean[3]), "&",
    f(btab$Married_mean[4]), "\\\\& (",
    f(btab$Married_sd[1]), ") & (",
    f(btab$Married_sd[2]), ") & (",
    f(btab$Married_sd[3]), ") & (",
    f(btab$Married_sd[4]), ") \\\\ \\\\
      Percent No College & ",
    f(btab$HS_mean[1]), "&",
    f(btab$HS_mean[2]), "&",
    f(btab$HS_mean[3]), "&",
    f(btab$HS_mean[4]), "\\\\& (",
    f(btab$HS_sd[1]), ") & (",
    f(btab$HS_sd[2]), ") & (",
    f(btab$HS_sd[3]), ") & (",
    f(btab$HS_sd[4]), ") \\\\
      \\bottomrule \\\\
      Number of Observations & \\multicolumn{1}{c}{",
    f(btab$n[1]), "} & \\multicolumn{1}{c}{",
    f(btab$n[2]), "} & \\multicolumn{1}{c}{",
    f(btab$n[3]), "} & \\multicolumn{1}{c}{",
    f(btab$n[4]), "}\\\\",
    sep = "",
    file = paste0(opath, "Tables/bal_table.tex")
)




# Financial Literacy ------------------------------------------------------


# subset dat to just group 1 and 2
dat12 <-
    dat %>%
    filter(group_assigned %in% c(1, 2))








# Data Manipulation -------------------------------------------------------
dat$risk_aversion <- ifelse(grepl("200", dat$la7), "Less Risk Averse", "More Risk Averse")
summary(dat$risk_aversion)

dat$loss_aversion <-
    ifelse(
        (
            (
                grepl("sure", dat$aq10) &
                    grepl("probability", dat$aq11)
            ) |
                (
                    grepl("sure", dat$aq12) &
                        grepl("probability", dat$aq13)
                )
        ),
        "More Loss Averse",
        "Less Loss Averse"
    )

table(dat$loss_aversion)

dat$time <- dat$finish - dat$start

dat$marg_rate_1 <- ifelse(is.na(dat$q1_c1), dat$q1_c2, dat$q1_c1)
summary(dat$marg_rate_1)

dat$marg_rate_2 <- ifelse(is.na(dat$q2_c1), dat$q2_c2, dat$q2_c1)
summary(dat$marg_rate_2)

dat$marg_rate_3 <- ifelse(is.na(dat$q3_c1), dat$q3_c2, dat$q3_c1)
summary(dat$marg_rate_3)

dat$marg_rate_4 <- ifelse(is.na(dat$q4_c1), dat$q4_c2, dat$q4_c1)
summary(dat$marg_rate_4)

# Binary Versions
dat$retire_time_bin <- ifelse((dat$age_retire - dat$age_current) > median(dat$age_retire - dat$age_current),
    "Near Retirement",
    "Not Near Retirement"
)
dat$retire_size_bin <- ifelse((dat$retire_plan) > median(dat$retire_plan),
    "Large Current Plan",
    "Small Current Plan"
)



summary(dat$change1)
summary(dat$change2)
summary(dat$change3)
summary(dat$change4)

.16 * (log(.4) - log(.2)) / .3



dat$gross_expense_ratio_random <- case_when(
    dat$gross_expense_ratio_random == 2 ~ .003,
    dat$gross_expense_ratio_random == 3 ~ .005,
    dat$gross_expense_ratio_random == 4 ~ .007
)

a <- lm(I(1 - dat$change1) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
b <- lm(I(1 - dat$change2) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
c <- lm(I(1 - dat$change3) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
d <- lm(I(1 - dat$change4) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()

cat("log(Expense Ratio):Base Group & ",
    f(a$coefficients[2, 1], 3), " & ",
    f(b$coefficients[2, 1], 3), " & ",
    f(c$coefficients[2, 1], 3), " & ",
    f(d$coefficients[2, 1], 3), " \\\\
  & (",
    f(a$coefficients[2, 3], 3), ") & (",
    f(b$coefficients[2, 3], 3), ") & (",
    f(c$coefficients[2, 3], 3), ") & (",
    f(d$coefficients[2, 3], 3), ") \\\\ \\\\
  log(Expense Ratio):Default Group& ",
    f(a$coefficients[3, 1], 3), " & ",
    f(b$coefficients[3, 1], 3), " & ",
    f(c$coefficients[3, 1], 3), " & ",
    f(d$coefficients[3, 1], 3), " \\\\
  & (",
    f(a$coefficients[3, 3], 3), ") & (",
    f(b$coefficients[3, 3], 3), ") & (",
    f(c$coefficients[3, 3], 3), ") & (",
    f(d$coefficients[3, 3], 3), ") \\\\ \\\\
  log(Expense Ratio):Distribution Group & ",
    f(a$coefficients[4, 1], 3), " & ",
    f(b$coefficients[4, 1], 3), " & ",
    f(c$coefficients[4, 1], 3), " & ",
    f(d$coefficients[4, 1], 3), " \\\\
  & (",
    f(a$coefficients[4, 3], 3), ") & (",
    f(b$coefficients[4, 3], 3), ") & (",
    f(c$coefficients[4, 3], 3), ") & (",
    f(d$coefficients[4, 3], 3), ") \\\\ \\\\
  log(Expense Ratio):Default and Distribution & ",
    f(a$coefficients[5, 1], 3), " & ",
    f(b$coefficients[5, 1], 3), " & ",
    f(c$coefficients[5, 1], 3), " & ",
    f(d$coefficients[5, 1], 3), " \\\\
  & (",
    f(a$coefficients[5, 3], 3), ") & (",
    f(b$coefficients[5, 3], 3), ") & (",
    f(c$coefficients[5, 3], 3), ") & (",
    f(d$coefficients[5, 3], 3), ") \\\\",
    file = paste0(opath, "tables/switching.tex")
)















l <- which(grepl("Dist", dat$group_assigned))

lm(change1 ~ I(log(gross_expense_ratio_random)) * loss_aversion + I(log(gross_expense_ratio_random)) * risk_aversion, data = dat) %>% summary()
lm(change2 ~ I(log(gross_expense_ratio_random)) * loss_aversion + I(log(gross_expense_ratio_random)) * risk_aversion, data = dat) %>% summary()
lm(change3 ~ I(log(gross_expense_ratio_random)) * loss_aversion + I(log(gross_expense_ratio_random)) * risk_aversion, data = dat) %>% summary()
lm(change4 ~ I(log(gross_expense_ratio_random)) * loss_aversion + I(log(gross_expense_ratio_random)) * risk_aversion, data = dat) %>% summary()


lm(I(q4_a == "Buffer RILA") ~ group_assigned * loss_aversion, data = dat) %>% summary()

# Initial Choice ----------------------------------------------------------
# Reshape the data for combined plot with fin_lit
dat_long <- dat %>%
    select(q1_a, q2_a, q3_a, fin_lit) %>%
    pivot_longer(
        cols = c(q1_a, q2_a, q3_a),
        names_to = "question_type",
        values_to = "value"
    )

# Ensure fin_lit is numeric
dat_long$fin_lit <- as.numeric(dat_long$fin_lit)

# Calculate percentages and confidence intervals with fin_lit grouping
combined_data <- dat_long %>%
    group_by(question_type, fin_lit, value) %>%
    summarise(
        count = n(),
        .groups = "drop"
    ) %>%
    group_by(question_type, fin_lit) %>%
    mutate(
        total = sum(count),       # Total counts per 'question_type' and 'fin_lit'
        perc = count / total,     # Percentage within each group
        ci = calc_ci(perc, total) # Confidence interval
    ) %>%
    ungroup()

# Create the plot with error bars and fin_lit faceting
p <- ggplot(combined_data, aes(x = value, y = perc)) +
    geom_col() +
    geom_errorbar(
        aes(
            ymin = pmax(0, perc - ci),
            ymax = pmin(1, perc + ci)
        ),
        width = 0.2
    ) +
    facet_grid(
        fin_lit ~ question_type,
        labeller = as_labeller(c(
            q1_a = "TDF - Initial Choice",
            q2_a = "TD-RILA with Floor - Initial Choice",
            q3_a = "TD-RILA with Buffer - Initial Choice",
            `1` = "High Financial Literacy",
            `0` = "Low Financial Literacy"
        ))
    ) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Investment Choice - Risk Level") 

# Display or save the plot
if (sa == 1) {
    ggsave(file = paste0(opath, "Figures/combined_plot_by_finlit.jpg"), 
           plot = p, 
           width = width, 
           height = height)
} else {print(p)
}



# q1a
dat %>%
    # group_by(group_assigned) %>%
    # mutate(count = n()) %>%
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>%
    ggplot(aes(x = q1_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) +
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") +
    xlab("TDF - Initial Choice") +
    theme(
        text = element_text(size = t1), # Global text size
        axis.title = element_text(size = t2), # Axis title size
        axis.text = element_text(size = t2) # Axis text size
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q1a1.jpg"),
        width = width,
        height = height
    )
}

dat %>%
    # group_by(group_assigned) %>%
    # mutate(count = n()) %>%
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>%
    ggplot(aes(x = q2_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) +
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") +
    xlab("TD-RILA with Floor - Initial Choice") +
    theme(
        text = element_text(size = t1), # Global text size
        axis.title = element_text(size = t2), # Axis title size
        axis.text = element_text(size = t2) # Axis text size
    )


if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q2a1.jpg"),
        width = width,
        height = height
    )
}


dat %>%
    # group_by(group_assigned) %>%
    # mutate(count = n()) %>%
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>%
    ggplot(aes(x = q3_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) +
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") +
    xlab("TD-RILA with Buffer - Initial Choice") +
    theme(
        text = element_text(size = t1), # Global text size
        axis.title = element_text(size = t2), # Axis title size
        axis.text = element_text(size = t2) # Axis text size
    )


if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q3a1.jpg"),
        width = width,
        height = height
    )
}

dat %>%
    filter(group_assigned == "1. Base") %>%
    # mutate(count = n()) %>%
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) +
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") +
    xlab("TDF - Initial Choice") +
    theme(
        text = element_text(size = t1), # Global text size
        axis.title = element_text(size = t2), # Axis title size
        axis.text = element_text(size = t2) # Axis text size
    )


if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a1.jpg"),
        width = width,
        height = height
    )
}

# Modified plot for q4a with confidence intervals
q4a_data <- dat %>%
    group_by(group_assigned) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, q4_a) %>%
    summarise(
        perc = n() / first(count),
        count = first(count),
        .groups = "drop"
    ) %>%
    mutate(ci = calc_ci(perc, count))


# Create the plot
ggplot(q4a_data, aes(x = q4_a, y = perc)) +
    geom_col() +
    geom_errorbar(aes(ymin = pmax(0, perc - ci), ymax = pmin(1, perc + ci)), width = 0.2) +
    facet_grid(. ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice") +
    theme(
        text = element_text(size = t1),
        axis.title = element_text(size = t2),
        axis.text = element_text(size = t2)
    )

# Save the plot if sa == 1
if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a.jpg"),
        width = width2,
        height = height
    )
}

# Modified plot for q4a by risk aversion with confidence intervals
q4a_ra_data <- dat %>%
    group_by(group_assigned, risk_aversion) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, risk_aversion, q4_a) %>%
    summarise(
        perc = n() / first(count),
        count = first(count),
        .groups = "drop"
    ) %>%
    mutate(ci = calc_ci(perc, count))

ggplot(q4a_ra_data, aes(x = q4_a, y = perc)) +
    geom_col() +
    geom_errorbar(aes(ymin = pmax(0, perc - ci), ymax = pmin(1, perc + ci)), width = 0.2) +
    facet_grid(risk_aversion ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice") +
    theme(
        text = element_text(size = t1),
        axis.title = element_text(size = t2),
        axis.text = element_text(size = t2)
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a_ra_with_ci.jpg"),
        width = width,
        height = height
    )
}



# q2a
dat %>%
    group_by(group_assigned) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, q2_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q2_a, y = perc)) +
    geom_col() +
    facet_grid(. ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Floor RILA - Initial Choice") +
    theme(
        text = element_text(size = t1), # Global text size
        axis.title = element_text(size = t2), # Axis title size
        axis.text = element_text(size = t2) # Axis text size
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q2a.jpg"),
        width = width2,
        height = height
    )
}

# q4a - ra
dat %>%
    group_by(group_assigned, risk_aversion) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, risk_aversion, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(risk_aversion ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice")

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a_ra.jpg"),
        width = width,
        height = height
    )
}


# q4a - la
dat %>%
    group_by(group_assigned, loss_aversion) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, loss_aversion, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(loss_aversion ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice")

# q4a - la
dat %>%
    group_by(group_assigned, loss_aversion) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, loss_aversion, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(loss_aversion ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice")

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a_la.jpg"),
        width = width,
        height = height
    )
}

# q4a - la
dat %>%
    group_by(group_assigned, retire_time_bin) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, retire_time_bin, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(retire_time_bin ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice")

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a_rt.jpg"),
        width = width,
        height = height
    )
}


# q4a - la
dat %>%
    group_by(group_assigned, retire_size_bin) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, retire_size_bin, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(retire_size_bin ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice")

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a_rs.jpg"),
        width = width,
        height = height
    )
}



dat %>%
    group_by(group_assigned, fin_lit) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, fin_lit, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(fin_lit ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice")


# Income-based analysis
dat$income_group <- ifelse(dat$salary_current > median(dat$salary_current),
    "High Income", "Low Income"
)

# Plot for income groups
dat %>%
    group_by(group_assigned, income_group) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, income_group, q4_a) %>%
    summarise(perc = n() / count) %>%
    distinct() %>%
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(income_group ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Combo - Initial Choice") +
    theme(
        text = element_text(size = t1), # Global text size
        axis.title = element_text(size = t2), # Axis title size
        axis.text = element_text(size = t2) # Axis text size
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q4a_income.jpg"),
        width = width,
        height = height
    )
}

# Modified plot for q1a with confidence intervals
q1a_data <- dat %>%
    group_by(group_assigned) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, q1_a) %>%
    summarise(
        perc = n() / first(count),
        count = first(count),
        .groups = "drop"
    ) %>%
    mutate(ci = calc_ci(perc, count))

ggplot(q1a_data, aes(x = q1_a, y = perc)) +
    geom_col() +
    geom_errorbar(aes(ymin = pmax(0, perc - ci), ymax = pmin(1, perc + ci)), width = 0.2) +
    facet_grid(. ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("TDF - Initial Choice") +
    theme(
        text = element_text(size = t1),
        axis.title = element_text(size = t2),
        axis.text = element_text(size = t2)
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q1a.jpg"),
        width = width2,
        height = height
    )
}

# Modified plot for q2a with confidence intervals
q2a_data <- dat %>%
    group_by(group_assigned) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, q2_a) %>%
    summarise(
        perc = n() / first(count),
        count = first(count),
        .groups = "drop"
    ) %>%
    mutate(ci = calc_ci(perc, count))

ggplot(q2a_data, aes(x = q2_a, y = perc)) +
    geom_col() +
    geom_errorbar(aes(ymin = pmax(0, perc - ci), ymax = pmin(1, perc + ci)), width = 0.2) +
    facet_grid(. ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Floor RILA - Initial Choice") +
    theme(
        text = element_text(size = t1),
        axis.title = element_text(size = t2),
        axis.text = element_text(size = t2)
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q2a.jpg"),
        width = width2,
        height = height
    )
}

# Modified plot for q3a with confidence intervals
q3a_data <- dat %>%
    group_by(group_assigned) %>%
    mutate(count = n()) %>%
    group_by(group_assigned, q3_a) %>%
    summarise(
        perc = n() / first(count),
        count = first(count),
        .groups = "drop"
    ) %>%
    mutate(ci = calc_ci(perc, count))

ggplot(q3a_data, aes(x = q3_a, y = perc)) +
    geom_col() +
    geom_errorbar(aes(ymin = pmax(0, perc - ci), ymax = pmin(1, perc + ci)), width = 0.2) +
    facet_grid(. ~ group_assigned) +
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") +
    xlab("Buffer RILA - Initial Choice") +
    theme(
        text = element_text(size = t1),
        axis.title = element_text(size = t2),
        axis.text = element_text(size = t2)
    )

if (sa == 1) {
    ggsave(
        file = paste0(opath, "Figures/q3a.jpg"),
        width = width2,
        height = height
    )
}

